{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <center> Открытый курс по машинному обучению. Сессия № 3 </center>\n",
    "#### <center> Индивидуальный проект по анализу данных </center>\n",
    "### <center> Предсказание дисфункции щитовидной железы </center>\n",
    "<center> Автор: Мирский Борис Юрьевич (boris.y.mirsky@gmail.com) </center> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Описание набора данных и признаков\n",
    "Для анализа был взят датасет http://archive.ics.uci.edu/ml/datasets/Thyroid+Disease с данными обследования пациентов по выявлению дисфункции щитовидной железы (далее используется сокращение '__щ.ж.__'). Этот набор данных датирован 1987 годом и был собран в Австралии. В развитых странах примерно 1% населения имеет проблемы с щ.ж., которые понижают качество жизни и способны перерасти в более серьёзные заболевания.\n",
    "\n",
    "Полностью датасет состоит из 6 наборов данных одного размера. Каждый набор заранее разбит на train и test, также есть отдельный файл с названиями признаков. Есть предупреждение о повреждении некоторых наборов и о пересечении части данных. Данное исследование ограничивается одним из шести доступных наборов данных. \n",
    "\n",
    "1. __age__  (численный):\t\t\t\t Возраст в годах. Возможные значения - float.\n",
    "2. __sex___ (бинарный):\t\t\t\tГол. Возможные значения - F,M\n",
    "3. __on thyroxine__ (бинарный):  \t ПРЕДПОЛОЖИТЕЛЬНО: принимает ли пациент синтетический аналог тироксина (тироксин - гормон, выделяемый в кровь щ. ж.). Возможные значения - t, f.\n",
    "4. __query on thyroxine__ (бинарный):\t\t Есть ли отклонения по уровню тироксина. Возможные значения - t, f.\n",
    "5. __on antithyroid medication__ (бинарный):\tУпотребление пациентом антихироидных препаратов. Возможные значения - t, f.\n",
    "6. __sick__ (бинарный):\t\t\t\tБолеет ли пациент (видимо, болеет чем-то помимо щ.ж.). Возможные значения - t, f.\n",
    "7. __pregnant__ (бинарный):\t\t\tБеременность. Возможные значения - t, f.\n",
    "8. __thyroid surgery__ (бинарный):\tБыло ли хирургическое вмешательство по поводу щ. ж. Возможные значения - t, f. \n",
    "9. __I131 treatment__ (бинарный):\tПроводилось ли лечение йодом I-131 (радиоактивный изотоп йода). Возможные значения - t, f.\n",
    "10. __query hypothyroid__ (бинарный): Есть ли у пациента состояние, вызванное нехваткой гормона, произодимого щ.ж. Возможные значения - t, f.\n",
    "11. __query hyperthyroid__ (бинарный):\tЕсть ли у пациента перепроизводство щ.ж. своего гормона.\n",
    "12. __lithium__ (бинарный):\t\t\t    Принимает ли пациент литий. Возможные значения - t, f.\n",
    "13. __goitre__ (бинарный):\t\t\t\tЕсть ли у пациента зоб (разрастание щ.ж.). Возможные значения - t, f.\n",
    "14. __tumor__ (бинарный):\t\t\t\tЕсть ли у пациента опухоль. Возможные значения - t, f.\n",
    "15. __hypopituitary__ (бинарный):\t\tЕсть ли у пациента снижение секреции одного или нескольких из восьми гормонов. Возможные значения - t, f.\n",
    "16. __psych__ (бинарный):\t\t\t\tПрибегал ли пациент к помощи психотерапевта. Возможные значения - t, f.\n",
    "17. __TSH measured__ (бинарный):\t\tИзмерялся ли Thyroid stimulating hormone (Tиреотропный гормон). Возможные значения - t, f.\n",
    "18. __TSH__ (численный):                Измеренный тиреотропный гормон. Возможные значения - float.\n",
    "19. __T3  measured__ (бинарный):\t\tИзмерялся ли гормон  Triiodothyronine (трийодтиронин). Возможные значения - t, f.\n",
    "20. __T3__ (численный):\t\t\t\t    Измеренный гормон трийодтирпациентонин. Возможные значения - float.\n",
    "21. __TT4 measured__ (бинарный):\t\tИзмерялся ли общий сывороточный тироксин (Total serum thyroxine). Возможные значения - t, f.\n",
    "22. __TT4__ (численный):\t\t\t\tИзмеренный гормон тироксин. Возможные значения - float. \n",
    "23. __T4U measured__ (бинарный):\t\tИзмерялся ли гормон\tсвободный тироксин (thyroxine). Возможные значения. - t, f.\n",
    "24. __T4U__ (численный):\t\t\t\tИзмеренный гормон. Возможные значения - float.\n",
    "25. __FTI measured__ (бинарный):\t\tИзмерялся ли индекс свободного тироксина ( Free Tyroxine Index). Возможные значения - t, f.\n",
    "26. __FTI__ (численный):\t\t\t\tИзмеренный индекса свободного тироксина. Возможные значения - float. \n",
    "27. __TBG measured__ (бинарный):\t\tИзмерялся ли  Thyroid binding globulin (тироксинсвязывающий глобулин). Возможные значения - t, f.      \n",
    "28. __TBG__ (численный):\t            Измеренный тироксинсвязывающий глобулин. Возможные значения - float.\n",
    "29. __referral source__ (категориальный):\tРекомендации: WEST (установить не удалось); STMW (пить меньше воды, если пить - только воду с минералами); SVHC (обязательные регулярные осмотры на thyroxin); SVI (если уровень тироксина возрастёт - срочно к врачу); SVHD (Требуется антитиреоидная (Antithyroid) терапия); others (другие). Возможные значения - строковые. \n",
    "30. __pred__ (категориальный):          __Целевая переменная__. Три состояния: increased binding protein, decreased binding proteinHyper, negative. Возможные значения - строковые. \n",
    "\n",
    "#### Цель данного проекта - предсказать дисфункцию щ.ж. у недиагностированных пациентов через вынесение одного из трёх диагнозов:\n",
    "1. increased binding protein: Hypo Thyroid (гипотериоз т.е. недостаток тероида) \n",
    "2. decreased binding proteinHyper: Hyper Thyroid (гипертериоз т.е. переизбыток тероида) \n",
    "3. negative: Normal (пациент здоров) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import random\n",
    "import itertools\n",
    "from matplotlib import pyplot as plt\n",
    "import seaborn as sns \n",
    "import pylab as plt\n",
    "from sklearn.utils import shuffle\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.ensemble import RandomForestClassifier, GradientBoostingClassifier, ExtraTreesClassifier\n",
    "from sklearn.model_selection import train_test_split, learning_curve, validation_curve\n",
    "from sklearn.model_selection import StratifiedShuffleSplit, GridSearchCV\n",
    "from sklearn.metrics import auc, roc_auc_score, roc_curve, precision_recall_curve, classification_report\n",
    "from sklearn.metrics import accuracy_score, confusion_matrix, f1_score\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "pd.options.mode.chained_assignment = None\n",
    "RANDOM_STATE=42\n",
    "%matplotlib inline\n",
    "sns.set_context(\"notebook\")\n",
    "plt.style.use('ggplot')\n",
    "plt.rcParams['figure.figsize']=(8, 6)  \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Первичный анализ данных \n",
    "## 3. Первичный визуальный анализ данных \n",
    "## 4. Инсайты, найденные зависимости\n",
    "## 7. Предобработка данных\n",
    "\n",
    "Дисфункции щ.ж. относятся к области медицины 'эндокринология', изучающей и лечащей гормональную систему человека. Поэтому часть признаков обозначает гормоны: данные анализов представлены аббревиатурами вида 'TT4'; признаки вида 'measured_TT4' говорят о самом факте анализа. \n",
    "\n",
    "Также есть обычные параметры: возраст; пол; беременность; было ли хирургическое вмешательств по поводу заболевания щ.ж.; есть ли опухоль; есть ли зоб (разросшаяся опухоль); наблюдается ли пациент у психотерапевта; принимает ли препараты лития (успокоительные) и т.д. \n",
    "\n",
    "Наш датасет уже разбит на train и test, что неудобно для некоторых действий по анализу и предобработке. Поэтому объединим обе части в одну выборку. Затем, после нескольких операций, снова разобьём на train и test."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_ = pd.read_csv('allbp.data.txt', header=None)                                   # (2800, 30)\n",
    "# names_ = pd.read_csv('allbp.names.txt', error_bad_lines=False, header=None)         # (13, 1)\n",
    "test_ = pd.read_csv('allbp.test.txt', header=None)                                    # (972, 30) \n",
    "\n",
    "# Это наш датасет из необработанных данных\n",
    "X_ = pd.concat([train_, test_])  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Файл names содержит названия признаков, но из-за особенностей форматирования оказался бесполезен. Руками извлекаем имена столбцов и присваиваем их нашей выборке."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "names = ['age','sex','on thyroxine','query on thyroxine','on antithyroid medication','sick','pregnant',\n",
    "      'thyroid surgery','I131 treatment','query hypothyroid','query hyperthyroid','lithium','goitre',\n",
    "      'tumor','hypopituitary','psych','TSH measured','TSH','T3 measured','T3','TT4 measured','TT4',\n",
    "      'T4U measured','T4U','FTI measured','FTI','TBG measured','TBG','referral source', 'pred']\n",
    "\n",
    "X_.columns = names "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hаши сырые данные"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_.head() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_.shape "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_.info(verbose=True, null_counts=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Функция info() показывает: (1) отсутствие пустых строк т.к. вместо NaN'ов стоят '?'; (2) принадлежность всех объектов к типу object. На самом деле признаки принадлежат двум типам: \n",
    "* Вещественные числа 'num_features' - 6 столбцов.\n",
    "* Категориальные признаки 'nom_features' - 24 столбца. Из них 22 принимают 2 значения, т.е. бинарные. \n",
    "* Столбец 'pred' это наша целевая переменная. \n",
    "\n",
    "Посмотрим, в каких столбцах есть строки с пропущенными данными. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nan_list = []\n",
    "for i in X_:\n",
    "    if len(X_[X_[i] == '?']) != 0:\n",
    "        nan_list.append((i, (len(X_[X_[i] == '?']))))\n",
    "print(nan_list) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Пропущенные значения распределены неравномерно: столбец TBG весь такой, у T3 их 21%, у остальных менее 10%. Oтметим, что список nan_list почти полностью (кроме 'sex') совпадает со списком num_features (ниже).\n",
    "\n",
    "#### Преобразование признаков - 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# создадим новый пустой DataFrame, в который будем заносить обработанные данные\n",
    "X = pd.DataFrame() \n",
    "\n",
    "# числовые признаки \n",
    "num_features = ['age', 'TSH', 'T3', 'TT4', 'T4U', 'FTI'] \n",
    "# категориальные признаки, все кроме pred\n",
    "# Только referral source настоящий номинальный - referral source, остальные бинарные. \n",
    "nom_features = ['sex', 'on thyroxine', 'query on thyroxine', 'on antithyroid medication', 'sick', 'pregnant', \n",
    "                   'thyroid surgery', 'I131 treatment', 'query hypothyroid', 'query hyperthyroid', \n",
    "                   'lithium', 'goitre', 'tumor', 'hypopituitary', 'psych', 'TSH measured', 'T3 measured',\n",
    "                   'TT4 measured', 'T4U measured', 'FTI measured', 'TBG measured', 'TBG', 'referral source'] \n",
    "\n",
    "# подготовим отображение \n",
    "nom_features_mapping = [{\"F\": 0, \"M\": 1},\n",
    "    {\"f\": 0, \"t\": 1}, {\"f\": 0, \"t\": 1}, {\"f\": 0, \"t\": 1}, {\"f\": 0, \"t\": 1}, {\"f\": 0, \"t\": 1},\n",
    "    {\"f\": 0, \"t\": 1}, {\"f\": 0, \"t\": 1}, {\"f\": 0, \"t\": 1}, {\"f\": 0, \"t\": 1}, {\"f\": 0, \"t\": 1},\n",
    "    {\"f\": 0, \"t\": 1}, {\"f\": 0, \"t\": 1}, {\"f\": 0, \"t\": 1}, {\"f\": 0, \"t\": 1}, {\"f\": 0, \"t\": 1},\n",
    "    {\"f\": 0, \"t\": 1}, {\"f\": 0, \"t\": 1}, {\"f\": 0, \"t\": 1}, {\"f\": 0, \"t\": 1}, {\"f\": 0, \"t\": 1},\n",
    "    {\"f\": 0, \"t\": 1}, {\"WEST\":0, \"STMW\":1, \"SVHC\":2, \"SVI\":3, \"SVHD\":4, \"other\":5}]  \n",
    "\n",
    "# Функция для замены '?' на NaN    \n",
    "nan_func = lambda x: np.nan if str(x) == \"?\" else x\n",
    "\n",
    "# Обработаем num_features\n",
    "for f in num_features:\n",
    "    X[f] = X_[f].map(nan_func).astype(\"float64\") \n",
    "     \n",
    "# Обработаем nоm_features, отобразив категориальные значения в численные\n",
    "for n, f in enumerate(nom_features):\n",
    "    X[f] = X_[f].map(nan_func).map(nom_features_mapping[n])  \n",
    "\n",
    "# Столбец 'pred' требует отдельной обработки.\n",
    "# Это категориальный признак с тремя классами, приведём его к численному виду.\n",
    "# Цифры после .| неважны, отбросим их. \n",
    "\n",
    "y_ = pd.Series([])\n",
    "for i in range(0, X.shape[0]):\n",
    "    y_[i] = X_['pred'].iloc[i].split('.')[0] \n",
    "# Отобразим категориальные значения на численные\n",
    "y = y_.map({\"negative\":0, \"increased binding protein\":1, \"decreased binding protein\":2})  \n",
    "X['pred'] = y \n",
    "\n",
    "# Заполнять пустые ячейки средними значениями будем позже, но для 'sex' заполним нулями сейчас   \n",
    "X.loc[:,'sex'].fillna(0, inplace=True) \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Распределение признаков\n",
    "Посмотрим на распределение целевой перемeнной "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.value_counts(X['pred'].values) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# То же самое более наглядно\n",
    "X['pred'].value_counts().plot(kind=\"bar\")     #, label=\"ckd\")\n",
    "plt.legend()\n",
    "plt.title(\"Distribution of classes\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Видим, что целевая переменная крайне разбалансирована: \n",
    "* negative 95.1% \n",
    "* increased binding protein 4.6%\n",
    "* decreased binding protein 0.3% "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Проверим распределение остальных категориальных признаков."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for feature in nom_features:\n",
    "    print(\"Feature: %s\" % feature, sorted(X[feature].value_counts())) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Целевая переменная и все категориальные признаки крайне расбалансированы. Вспомним предупреждение о низкой доле таких больных: для этой области медицины это нормальная ситуация, и с ней тоже надо уметь работать. Признак 'sex', единственный без огромного разрыва, всё равно имеет более чем двойной перекос. Это наблюдение подтверждает справочные данные - женщины гораздо чаще подвержены нарушениям эндокринной системы. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Преобразование признаков - 2 (Features selection)\n",
    "\n",
    "Выше мы заметили, что\n",
    "* в стоолбце 'TBG' все значения = NaN. \n",
    "* два признака - 'TBG measured' и 'hypopituitary' - распределены [1, 3771] и [0, 3772] соответственно. \n",
    "\n",
    "Посмотрим на пары вида 'TSH measured' и 'TSH'. Здравый смысл подсказывает, что 5 признаков, начинающихся со слова 'measured', говорят только о самом факте измерения одного из гормонов, не несут никакой полезной информации и будут сильно коррелировать со своими же данными по уровню гормонов. Если мы сделаем такую выборку и посмотрим хотя бы на 50 строк, то увидим полное совпадение между 0 в признаке вида 'measured' и NaN'ом в столбце со значением гормона. \n",
    "\n",
    "Проверим нашe предположение ещё одним способом. Сделаем матрицу корреляций, включив туда претендентов на удаление, показатели гормонов в качестве маркеров и целевую переменную."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_corr_ = X[['TBG', 'TBG measured', 'hypopituitary', 'TSH measured', 'TSH', 'T3 measured', 'T3', 'TT4 measured', 'TT4', 'T4U measured', 'T4U', 'FTI measured', 'FTI', 'pred']].corr()    \n",
    "sns.heatmap(X_corr_, annot=True, cmap=\"Blues\") "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Видим серую (т.е. пустую) полосу слева и сверху таблицы, а также пустые ячейки вдоль диагонали. Наша гипотеза подтверждается. \n",
    "\n",
    "Удалим 8 бесполезных признаков. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = X.drop(['TBG', 'TBG measured', 'hypopituitary', 'TSH measured', 'T3 measured', 'TT4 measured', 'T4U measured', 'FTI measured'], axis=1) \n",
    "X.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Разобьём на train/test.\n",
    "\n",
    "Т.к. в численных признаках пропуски заменяются средним значением, деление на train/test выполним сейчас, чтобы не использовать данные из test для подсчёта среднего по признаку. Разбиение будет стратифицированным - это крайне важно при работе с несбалансированными данными. В функции разбиения параметр shuffle=True по умолчанию. На тестовую часть выделим 25% датасета (изначально было 26.8%)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = X.drop('pred', 1)\n",
    "\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25, random_state=42, stratify = y) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Пропуски в X_train[num_features] и X_test[num_features] заполним средними значениями."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for feature in num_features:\n",
    "    mean_val = X_train[feature].mean()                    \n",
    "    X_train[feature].fillna(mean_val, inplace=True, axis=0) \n",
    "    \n",
    "for feature in num_features:\n",
    "    mean_val = X_test[feature].mean()                    \n",
    "    X_test[feature].fillna(mean_val, inplace=True, axis=0) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "После того, как вычистили NaN'ы, можно посмотреть на распределение численных признаков.\n",
    "В левом столбце значения признаков, в правом они же прологарифмированные."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nf_count = len(num_features)\n",
    "_, axes = plt.subplots(nf_count, 2, figsize=(12, nf_count * 4))\n",
    "\n",
    "for i, col in enumerate(num_features):\n",
    "    sns.distplot(X_train[col], ax=axes[i, 0]);  \n",
    "    sns.distplot(np.log(X_train[col] + 1), ax=axes[i, 1]);      "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Все численные признаки имеют нормальное распределение с небольшой дисперсией. Напомню, что мы заполнили средним  21% у Т3, и по 5-8% у TSH, TT4, T4U, FTI. Видим, что эта правка нигде не исказила данные. \n",
    "\n",
    "Единственный признак с явными выбросами это 'age', где примерно в 38, 58 и 78 лет дисфункция явно учащается. Это связано с климаксом (возрастной перестройкой организма), протекающим в несколько стадий. \n",
    "\n",
    "Судя по графику, надо прологарифмировать только один признак - TSH. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train['TSH_log'] = np.log(X_train['TSH'])\n",
    "X_train = X_train.drop('TSH', axis=1) \n",
    "\n",
    "X_test['TSH_log'] = np.log(X_test['TSH'])\n",
    "X_test = X_test.drop('TSH', axis=1) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Преобразования признаков - 3\n",
    "После заполнения NaN'ов мы может провести нормализацию численных признаков. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scaler = StandardScaler()\n",
    "\n",
    "X_train[['age', 'T3', 'TT4', 'T4U', 'FTI', 'TSH_log']] = scaler.fit_transform(X_train[['age', 'T3', 'TT4', 'T4U', 'FTI', 'TSH_log']])\n",
    "X_test[['age', 'T3', 'TT4', 'T4U', 'FTI', 'TSH_log']] = scaler.transform(X_test[['age', 'T3', 'TT4', 'T4U', 'FTI', 'TSH_log']]) \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Преобразуем категориальный признак 'referral source' с помощью dummy кодирования."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train = pd.get_dummies(X_train, columns=['referral source']) \n",
    "X_test = pd.get_dummies(X_test, columns=['referral source']) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Основную работу по  __feature selection__ мы провели выше:\n",
    "* Преобразование признаков - 1. Поменяли '?' на NaN. Отобразили численные значения на категориальные. Заменили пропуски в 'sex'.\n",
    "* Преобразование признаков - 2. Удалили 8 признаков. Разбили на train/test. Заменили NaN'ы на среднее у категориальных признаков. Прологарифмировали один численный признак.\n",
    "* Преобразование признаков - 3. Провели нормализацию численных признаков. Для 'referral source' сделали dummy кодирование.\n",
    "\n",
    "Теперь надо снизить размерность признакового пространства. Признаков у нас осталось не очень много, поэтому использование мощных методов типа PCA будет избыточным. Мы можем обойтись измерением 'важности признаков'. \n",
    "\n",
    "Из расшифровки сокращений значений признака 'referral_source' понятно, что это набор рекомендаций. Неясно, когда они даются (после или до вынесения диагноза) - разница может быть оказаться важной. Но в любом случае не будем исключать эти 5 столбцов."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# При выполнии отборa признаков с помощью случайного леса рекомендуется небольшое количество деревьев.\n",
    "etc = ExtraTreesClassifier(n_estimators=50, random_state=0)\n",
    "etc.fit(X_train[X_train.columns].values, y_train.values) \n",
    "importances = etc.feature_importances_\n",
    "std = np.std([tree.feature_importances_ for tree in etc.estimators_], axis=0)\n",
    "indices = np.argsort(importances)[::-1]\n",
    "\n",
    "for f in range(X_train.shape[1]):\n",
    "    print(\"%d. feature %d (%f) %s\" % (f+1 , indices[f], importances[indices[f]] , X_train.columns[indices[f]]))   \n",
    "\n",
    "plt.figure()\n",
    "plt.bar(range(X_train.shape[1]), importances[indices], color=\"r\", yerr=std[indices], align=\"center\") \n",
    "plt.xticks(range(X_train.shape[1]), indices)\n",
    "plt.xlim([-1, X_train.shape[1]])\n",
    "plt.show() "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Между 7-м и 8-м признаками видим 2-х кратную разность, пусть она будет естественной границей. Оставим только 7 первых признаков: T4U, T3, TT4, FTI, pregnant, TSH_log, age.\n",
    "\n",
    "Итак, самыми важными оказались показания гормонов, беременность и возраст. Из них шесть численные и один бинарный. \n",
    "\n",
    "__C этой выборкой мы будем работать дальше__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train_final = X_train[['T4U', 'T3', 'TT4', 'FTI', 'pregnant', 'TSH_log', 'age']]\n",
    "X_test_final = X_test[['T4U', 'T3', 'TT4', 'FTI', 'pregnant', 'TSH_log', 'age']]   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train_final.shape, y_train.shape, X_test_final.shape, y_test.shape "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 9. Создание новых признаков \n",
    "Это сложная узкоспециализировання область знаний. Медицинских знаний либо стороннего эксперта у меня нет, а интуиция ничего не подсказала. Единственное, что пришло в голову - 'поздняя беременность', которая осложняется с возрастом и, возможно, влияет на целевую переменную. Но таких женщин оказалось только 9. В итоге я понял, что не могу корректно создать новый признак, поэтому этап Feature Engenireeng пропущен. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Выбор метрики\n",
    "Мы не можем пропускать (false negative) больных людей, которых и так крайне мало (оба миноритарных класса), поэтому нас больше всего интересует полнота (recall). Это мнение высказывает, например, К.Воронцов. Специалисты, больше настроенные на решение практических задач (Ю.Кашницкий и А.Дьяконов), говорят о f1-мере и ROC AUC. Через функцию sklearn.metrics.classification_report мы будем получать все названные метрики, отдавая предпочтение гармоническому среднему f1. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Выбор модели "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Перед нами стоит медицинская задача, требующая однозначной интерпретации решений, поэтому сразу напрашивается термин 'деревья'.\n",
    "\n",
    "В случае с несбалансированными данными можно использовать\n",
    "* Библиотеку imbalanced learn (методы undersampling, oversampling и смешанные сложные методы типа SMOTE)\n",
    "* Линейные модели в комбинации (pipeline) с деревьями или бустингом.\n",
    "* Первые 2 возможности мы проигнориуем из-за ограничений на размер проекта. Остановимся только на LogisticRegression, RandomForrest и GradientBoostingClassifier."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Построим предсказания на отобранных классификаторах без настройки параметров."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "forest = RandomForestClassifier(n_estimators=50, random_state=17, class_weight='balanced')\n",
    "logit = LogisticRegression(random_state=17, class_weight= 'balanced', multi_class='ovr')   \n",
    "boosting = GradientBoostingClassifier(random_state=17) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "forest.fit(X_train_final, y_train)\n",
    "print('RandomForestClassifier')\n",
    "print(classification_report(y_test, forest.predict(X_test_final))) \n",
    "logit.fit(X_train_final, y_train)\n",
    "print('LogisticRegression')\n",
    "print(classification_report(y_test, logit.predict(X_test_final))) \n",
    "boosting.fit(X_train_final, y_train)\n",
    "print('GradientBoostingClassifier')\n",
    "print(classification_report(y_test, boosting.predict(X_test_final))) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "GradientBoostingClassifier ловит миноритарные классы лучше всех."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8. Кросс-валидация и настройка гиперпараметров модели"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Стратегия кросс-валидации и подбор параметров."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sss = StratifiedShuffleSplit(n_splits=5, test_size=0.3, random_state=17) \n",
    "\n",
    "logit_parameters = {'C': np.logspace(0, 1, 10)}\n",
    "\n",
    "boosting_parameters = {'n_estimators': [3, 5, 10, 15, 30],\n",
    "                       'learning_rate': [0.5, 0.8, 1],\n",
    "                       'max_depth': [3, 5, 7],\n",
    "                       'max_features': [0.3, 0.5, 1], \n",
    "                       'min_samples_leaf': [3, 5, 7],\n",
    "                                             }                                                      \n",
    "\n",
    "forest_parameters = {'n_estimators': [5, 10, 50],\n",
    "                     'max_depth': [3, 5, 8],\n",
    "                     'max_features': [0.5, 0.75, 1.0]} "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def learning_model(model, X_train, y_train, X_test, y_test):\n",
    "    model.fit(X_train, y_train)\n",
    "    y_pred = model.predict(X_test)\n",
    "    y_roc = model.predict_proba(X_test)[:,1]\n",
    "    score_f1 = f1_score(y_test, y_pred, average='weighted')\n",
    "    print('f1_score = {}'.format(score_f1))\n",
    "    return y_roc, score_f1 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid_logit = GridSearchCV(logit, logit_parameters, n_jobs=1, scoring ='f1_weighted', cv=sss)            \n",
    "grid_boosting = GridSearchCV(boosting, boosting_parameters, n_jobs=1, scoring ='f1_weighted', cv=sss)\n",
    "grid_forest = GridSearchCV(forest, forest_parameters, n_jobs=1, scoring ='f1_weighted', cv=sss)         "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_score = [] \n",
    "f1_scores = [] \n",
    "for grid in (grid_logit, grid_boosting, grid_forest):\n",
    "    y, f1 = learning_model(grid, X_train_final, y_train, X_test_final, y_test);\n",
    "    y_score.append(y)\n",
    "    f1_scores.append(f1) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Лучший результат показал Boosting, с ним и будем работать дальше. \n",
    "\n",
    "Посмотрим на параметры у нашего чемпиона"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid_boosting.best_params_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "При попытке обучить экземпляр класса с этими параметрами метрики ухудшились даже в сравнении с параметрами по умолчанию."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "boost_best = GradientBoostingClassifier(learning_rate=0.8, max_depth=5, max_features=0.3, \n",
    "                                        min_samples_leaf=3, n_estimators=15, random_state=17)\n",
    "boost_best.fit(X_train_final, y_train)\n",
    "print(classification_report(y_test, boost_best.predict(X_test_final))) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " Вручную удалось подобрать другие значения, неплохо улучшившие recall и f1-score. Этот алгоритм мы используем далее."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "boost_test = GradientBoostingClassifier(learning_rate=0.8, max_depth=3, max_features=1, \n",
    "                                        min_samples_leaf=3, n_estimators=3, random_state=17)\n",
    "boost_test.fit(X_train_final, y_train)\n",
    "print(classification_report(y_test, boost_test.predict(X_test_final))) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 10. Построение кривых валидации и обучения"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "max_feat = 1\n",
    "min_leaf = 3\n",
    "n_est = 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_learning_curve(estimator, title, X, y, label_y, ylim=None, cv=None, scoring='roc_auc', \n",
    "                        n_jobs=1, train_sizes=np.linspace(.1, 1.0, 10)):\n",
    "    plt.figure()\n",
    "    plt.title(title)\n",
    "    if ylim is not None:\n",
    "        plt.ylim(*ylim)\n",
    "    plt.xlabel(\"Training examples\")\n",
    "    plt.ylabel(label_y)\n",
    "    train_sizes, train_scores, test_scores = learning_curve(\n",
    "        estimator, X, y, cv=cv, n_jobs=n_jobs, train_sizes=train_sizes, scoring=scoring)\n",
    "    train_scores_mean = np.mean(train_scores, axis=1)\n",
    "    train_scores_std = np.std(train_scores, axis=1)\n",
    "    test_scores_mean = np.mean(test_scores, axis=1)\n",
    "    test_scores_std = np.std(test_scores, axis=1)\n",
    "    plt.grid()\n",
    "\n",
    "    plt.fill_between(train_sizes, train_scores_mean - train_scores_std,\n",
    "                     train_scores_mean + train_scores_std, alpha=0.1,\n",
    "                     color=\"r\")\n",
    "    plt.fill_between(train_sizes, test_scores_mean - test_scores_std,\n",
    "                     test_scores_mean + test_scores_std, alpha=0.1, color=\"g\")\n",
    "    plt.plot(train_sizes, train_scores_mean, 'o-', color=\"r\",\n",
    "             label=\"Training score\")\n",
    "    plt.plot(train_sizes, test_scores_mean, 'o-', color=\"g\",\n",
    "             label=\"Cross-validation score\")\n",
    "\n",
    "    plt.legend(loc=\"best\")\n",
    "    return plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_learning_curve(boost_test, 'Gradient Boosting', X_train_final, y_train, scoring='f1_weighted', \n",
    "                    cv=sss, n_jobs=1, label_y='f1')      "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_validation_curve(estimator, X, y,  cv_param_name, cv_param_values,ylim=None, cv=None,\n",
    "                        scoring='roc_auc', n_jobs=-1, train_sizes=np.linspace(.1, 1.0, 10)):\n",
    "    plt.figure()\n",
    "    if ylim is not None:\n",
    "        plt.ylim(*ylim)\n",
    "    plt.xlabel(cv_param_name) \n",
    "    plt.ylabel(\"Score\")\n",
    "    train_sizes, train_scores, test_scores = learning_curve(\n",
    "        estimator, X, y, cv=cv, n_jobs=n_jobs, train_sizes=train_sizes, scoring=scoring)\n",
    "    val_train, val_test = validation_curve(estimator, X, y, cv_param_name,\n",
    "                                           cv_param_values, cv=cv, scoring=scoring)\n",
    "    val_train_mean = np.mean(val_train, axis=1)\n",
    "    val_train_std = np.std(val_train, axis=1)\n",
    "    val_test_mean = np.mean(val_test, axis=1)\n",
    "    val_test_std = np.std(val_test, axis=1)\n",
    "    plt.grid()\n",
    "\n",
    "    plt.fill_between(cv_param_values,  val_train_mean -  val_train_std,\n",
    "                      val_train_mean +  val_train_std, alpha=0.1, color=\"b\")\n",
    "    plt.fill_between(cv_param_values, val_test_mean - val_test_std,\n",
    "                     val_test_mean + val_test_std, alpha=0.1, color=\"y\")\n",
    "    plt.plot(cv_param_values, val_train_mean, 'o-', color=\"b\",\n",
    "             label=\"Training score\")\n",
    "    plt.plot(cv_param_values, val_test_mean, 'o-', color=\"y\",\n",
    "             label=\"Cross-validation score\")\n",
    "\n",
    "    plt.legend(loc=\"best\")\n",
    "    return plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Зависимость качества от количества деревьев."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "estimators = np.arange(1, 40, 5)                     \n",
    "plot_validation_curve(boost_test, X_train_final, y_train, cv_param_name='n_estimators', \n",
    "                      cv_param_values= estimators,  scoring='f1_weighted')  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 11. Прогноз для тестовой выборки\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_confusion_matrix(cm, classes, normalize=False, title='Confusion matrix', cmap=plt.cm.Blues):\n",
    "    plt.imshow(cm, interpolation='nearest', cmap=cmap)\n",
    "    plt.title(title)\n",
    "    plt.colorbar()\n",
    "    tick_marks = np.arange(len(classes))\n",
    "    plt.xticks(tick_marks, classes, rotation=45)\n",
    "    plt.yticks(tick_marks, classes)\n",
    "    if normalize:\n",
    "        cm = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis]\n",
    "        print(\"Normalized confusion matrix\")\n",
    "    else:\n",
    "        print('Confusion matrix, without normalization')\n",
    "    thresh = cm.max() / 2.\n",
    "    for i, j in itertools.product(range(cm.shape[0]), range(cm.shape[1])):\n",
    "        plt.text(j, i, cm[i, j],\n",
    "                 horizontalalignment=\"center\",\n",
    "                 color=\"white\" if cm[i, j] > thresh else \"black\")\n",
    "    plt.tight_layout()\n",
    "    plt.ylabel('True label')\n",
    "    plt.xlabel('Predicted label')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "font = {'size' : 15}\n",
    "plt.rc('font', **font)\n",
    "cnf_matrix = confusion_matrix(y_test, boost_test.predict(X_test_final))  # grid_boosting\n",
    "plt.figure(figsize=(8, 6))\n",
    "plot_confusion_matrix(cnf_matrix, classes=['Normal', 'Hyper Thyroid', ' Hypo Thyroid '], title='GradientBoosting confusion matrix')\n",
    "plt.show() "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 12. Выводы\n",
    "Построены модели предсказания нахождения пациентов с дисфункцией щ.ж. Модель градиентного бустинга оказалась лучшей из рассмотренных. Подбор параметров по решётке с задачей не справился, при ручном подборе был получен более сильный алгоритм (см. In [177]), который использовался дальше. Полученный результат по метрике recall 99%/57%/75% можно назвать средним, что объясняется сложностью задачи (напомню, что миноритарные классы составляют 4.3% и 0.3%).\n",
    "\n",
    "В качестве возможностей улучшения прогнозирования можно посоветовать:\n",
    "* Консультации медика для понимания имеющихся признаков и возможного создания новых.\n",
    "* Надо больше признаков, причём скорее всго можно обойтись без дополнительных исследований. Например, вытащить информацию из имеющихся общих анализов, наследственности, местности проживания, пищевых либо бытовых привычек и т.д.\n",
    "* По графику кривых обучения видно, что их графики сходятся и увеличение данных добавит качества.\n",
    "* Продолжить эксперименты с параметрами GradientBoostingClassifier. Надо, в частности, найти баланс между обнаружением двух миноритарных классов.\n",
    "* Использовать библиотеку imbalanced learn, особенно методы undersampling. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
